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A probabilistic analysis of the direct simulation of a homogeneous gas is given. A 
hierarchy of equations similar to the BBGKY hierarchy for the reduced probability 
densities is derived. By invoking the molecular chaos assumption, an equation sim- 
ilar to the Boltzmann equation for the single particle probability density and the 
corresponding H-theorem is derived. 

I. INTRODUCTION 

Direct simulation Monte Carlo method (DSMC) is a standard method for solving the 
Boltzmann equation numerically. In this method space is divided into cells of volume AV^ 
and a large number of "particles" (A^ = 10^-10®) represent the real gas molecules. The 
evolution of the gas for a short time At is calculated in two steps. In the first step all 
particles are propagated for a time At without collisions. In the second step some randomly 
chosen pairs of particles in the same cell are allowed to collide and change their velocities 
without changing their positions. Number of pairs {n) chosen to make collision attempts is 
given by the formula n = RN'^At/2V where i? is a parameter we choose, is the number of 
particles in the cell and V is the volume of the cell. We call n number of collision attempts 
because not every chosen pair makes a collision. A pair is allowed to make a collision with a 
probability uar/R, where u is their relative velocity and ctt is the total cross section. The 
results are not sensitive to value of the parameter R as long as it is big enough such that 
very few pairs violate the condition ucr/R < 1 since average number of successful collision 
attempts 

^ (uar) _ RN^At (uar) _ {uar) 
^ R ~ 2V R~ ~ 2V ' 
is independent of R. Here {uaT) is the average of uot over all possible pairs. Although 

taking a very big R is acceptable theoretically, for practical reasons R should not chosen be 

too big either. 

The original method is due mainly to G. A. Bird. A seminal paper- of Bird gave some 
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heuristic arguments to justify its use. There are many good references on the subject. Ref. 
y has a good tutorial on the subject and Ref. Is] is a monograph on the subject by Bird 
himself which is a complete reference for the developments up to its publication year 1994. 
Also books on rarified gas dynamics devote many chapters to the subject and Ref. 14 and 
Ref. la are useful references in this category. 

A variant of the method was derived by Nanbu^ starting from the Boltzmann equation. 
To represent the evolution of the real gas such methods should converge to the true solution 
of the Boltzmann equation in the limit — ^ oo, AV 0, and At 0. Convergence proofs 
were given by Babovsky^ and Babovsky and lUner- for Nanbu's method and by Wagner- for 
Bird's method. 

For the evolution of the velocity distribution of a spatially uniform gas there is no need 
to divide physical space into cells and we can just work in velocity space. Although Bird 
recommended^ dividing real space into cells for studying the spatially homogeneous gas, we 
will show that this division is unnecessary. If we consider velocity space only and collide 
random pairs, we should obtain the evolution of the velocity distribution. The purpose of 
this paper is to study this stochastic process. 

These efforts to solve the Boltzmann equation using stochastic methods were driven by 
scientific applications and there was no motivation to use them as a pedagogical tool. It 
is surprising that similar stochastic algorithms for the homogeneous gas were conceived by 
people interested in using them as a pedagogical tool to demonstrate the evolution of a gas 
to the Maxwell-Boltzmann distribution. The earliest of such articles of which the author is 
aware is that of Novak-Bortz^° who studied the evolution of a gas of two-dimensional disks. 
Their algorithm is based on taking random pairs and colliding them with a probability pro- 
portional to U(7. Eger and Kressii modified this algorithm and Bonomo and Riggil2. applied 
the modification to hard disks. There are also other papers^^ii^ that do not use DSMC 
type stochastic processes to demonstrate the Maxwell-Boltzmann distribution. Although 
the DSMC method was well known, these papers do not reference papers on DSMC. Appar- 
ently the idea of stochastic methods for the evolution of a homogeneous gas was conceived 
for pedagogical applications independently. 

As mentioned, the DSMC algorithm can be used to demonstrate the approach of a velocity 
distribution to the Maxwell-Boltzmann distribution. The algorithm also gives an estimate 
of how many collisions is required to reach equilibrium and how various parameters affect 
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the evolution of the system. 

Although direct simulation is intuitively appealing, it is not clear that direct simulation 
algorithms represent the evolution of a real gas. The convergence proofs we have cited are 
formal and difficult to read. In this paper we prove that in the direct simulation appropriately 
normalized single particle probability distribution satisfies the Boltzmann equation for a 
homogeneous gas. The proof is relatively easy and intuitively appeahng and also its language 
is familiar to the physicist from the well known BBGKY hierarchy.— 

In Sec. II we consider the stochastic algorithm for a homogeneous gas. We derive a 
hierarchy of equations for the probability distribution of particles similar to the BBGKY 
hierarchy.^ We use the molecular chaos assumption to derive an equation similar to the 
Boltzmann equation for the single particle probability distribution /(v). We derive an H- 
theorem for /(v) and prove convergence to equilibrium. We also show how the equation for 
/(v) reduces to the Boltzmann equation for a particular choice of collision probabilities and 
derive Bird's "time counter" and "no time counter" methods. 



II. ANALYSIS OF THE DIRECT SIMULATION ALGORITHM FOR A 

HOMOGENEOUS GAS 

Consider a homogeneous gas of ^ 1 molecules without internal degrees of freedom. We 
randomly select pairs of molecules to collide. All possible pairs have an equal probability 
of 2/{N — 1)N to be selected. Suppose the velocities of the pair are and v^. The 
conditional probability that after the collision they have the velocities and v^i in the 
intervals (P^c and d?^D is T{yA, v^; v^;, v/?) (P^c d^^c- (From now on we will denote (Pw 
as dw for simplicity.) We also assume the symmetries 

T{wa, vb; vc, vd) = T(vc, ^d] v^, v^) (la) 
T{^A, vb; Vc, v^) = T(vb, v^; v^, vc). (lb) 

The total probability is unity and therefore 

j T{va, vb; Vc, v^,) dvc d\D = J T{\a, v^; vc, v^,) d\A d^B = 1- (2) 

Every selected pair makes a collision, although as we will show, by defining T{va, v^; vc, v^j) 
some of the collisions do not change velocities. After each collision new velocities of the 
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molecules are replaced by the old ones, and we select a new pair for the next collision. Of 
course there is the possibility of choosing the same pair with a very small probability. If that 
happens we let them collide again. We don't keep record of pairs that have made collisions 
already. 

We define f^^^vi, V2, . . . , vtv) as the probability density for the molecules. Because the 
molecules are indistinguishable, we require that the f^^^ be totally symmetric: 

/(^)(vi, . . . , Vi, . . . , Vj, . . . , Vjv) = /^^^Vl, . . . , Vj, . . . , Vi, . . . , vtv). (3) 

We also define the reduced probabihty densities 

/W(vi, V2, . . . , Va/) = j /^^^vi, V2, . . . , Vat) dwM+id^M+2 ■ ■ ■ dw m ■ (4) 
Because we will be dealing with pairs of particles, it is useful to define 

fi^P{^A,^B) = f^'^^H^l, . . . , Vi = VA, . . . , Vj = Vb, . . . , Vm). (5) 

That is, the velocities of the i,j pair are replaced by \^a,^b in the /(^'^^(vi, V2, . . . , vjv/) 
where i,j < M. We will also use the notation f^^^^-v; n) for /(*^)(vi, V2, . . . , vm) after the 
n^^ collision. 

The function f^^\\]n) satisfies the equation 

1 ^ ^ f N 

f^\v;n + l) = f^Ti^A,VB■,n)T{^rA,VB■,Vi,^rJ)dvAd^rB. (6) 

The meaning of Eq. is clear. If i,j is the last pair of molecules that has collided, then 
the probability of having Vj, Vj pairs after the collision is the probability of having initial 
velocities v^, (represented by f-j\vA, v^)) multiplied by the probability of ending with 
Vj,Vj (represented by T(va, v^; Vj, Vj)). The sum over i,j and the factor 1/N{N — 1) in 
Eq. ([6]) represents the fact that all pairs are possible with probability 1/N{N — 1). 

If we integrate Eq. (j6]) over \m+i,^m+2, ■ ■ ■ ,^n, we obtain 
/<«)(v;„+l) = (iV-M)(iV-M-l) 

^2{N-M) 



37Yr ^ / 4m+1^(va, vs; n)T{\A, vs; v^, vm+i) dvAdvs d\M+i 
' i=i ^ 



AI M 



+ iV(iV- 1)^ 4f^(vA,VB;n)T(vA,VB;Vi,Vj)dvArfvB. (7) 



5 



The f'^^\\]n + 1) depends on /(*^+^^(v; n); Eq. ^ represents a hierarchy of equations 
similar to the BBGKY hierarchy.— 
The first equation in the hierarchy is 

/«(v;n + l) = (l-2/iV)/«(v;n) 

+ ^y" /^^^Vyi, VB;n)T(vA, vb; vc, v)rfv^dvijdvc. (8) 

If we make the assumption of molecular chaos 

/(2)(v^,VB;n) = /«(v^;n)/«(vB;n), (9) 

we obtain a nonlinear equation for f^-^\v;n) similar to the Boltzmann equation. For large 
this approximation is almost exact as shown by the following argument. The velocities 

vi, V2 can be correlated only if particles one and two have collided with each recently. But 

this probability is of order which implies that for large A^, the velocity distributions of 

any two particles are uncorrelated. Present personal computers can handle = 10^-10^ so 

the assumption is almost exact. 

The key assumption in the argument for the validity of Eq. ([9]) is "recently." Two particles 

might be correlated for a short time, but after they have made a few collisions with other 

particles the correlations are expected to disappear. 

Another simplification occurs for large A^. The factor of 2/A^ in Eq. ([8]) is small and thus 

we can consider r = 2r;,/A^ to be a continuous parameter which we call the collision time. 

Then At = 2/N and [f^^\v;n + l) - /(^^(v; n)] /Ar can be written as df^^\v;T)/dT and 

Eq. (IE]) becomes 

^i^^^S^ = -/(i)(v;r) + 1 /«(v^;r)/«K;r)r(v^,v^; vc,v)rfv^t/v5rfvc. (10) 

From now on we will suppress the superscript (1) and the collision time r in f^^\v,T). 

Equation (|TOl) can be expressed as 

df(v) f 

= -fiy) + / fiyA)fiyB)T{yA, v^; vc, v) c?va d^B dvc- (11) 

A. The H-theorem and approach to equilibrium 

By using the relation 

/(v) = j f{v)f{^rc)T{vA,VB;vc,^)dvAdvBdvc, (12) 
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which follows from Eq. ([2]) and the normalization of f{"Vc), we can write Eq. (fTT!) as 

= j [/(va)/(vb) - /(v)/(vc)] T(va, Vb; Vc, v) dv^ rfv^ rfvc- (13) 

This form is similar to the Boltzmann equation. 

We can derive an H-theorem for this equation. We define H[t) as 

^(r)=y"/(v)ln(/(v))dv, (14) 

and use Eqs. ([1]) and (fT3|) to express dH/dr as 

^ = ~J -<il[f]T{vA,VB;vc,v)dvAdvBdvcdv, (15) 

where 

W] = [/(va)/(vb) - /(v)/(vc)] [ln/(v^)/(vB) - In /(v)/(vc)] . (16) 
The function can be shown to be always nonnegative. We argue that [x — y){\nx — Iny) 
is nonnegative for all positive x and y; Inx is an increasing function and thus x — y and 
In a; — In ?/ always have the same sign. Their product is always either positive or zero and 
zero occurs for x = y. T(v^, v^; vc, v) is intrinsically positive. Therefore the integrand is 
positive and dH/ dr is negative. 

Following the usual arguments of the H-theorem, the decrease of H stops only when 

In /(va) + In /(vs) = In /(v) + In /(vc) (17) 

is satisfied, which implies that ln/(v) is a collision invariant. If we choose T(va, v^; vcv) 
such that the total momentum and energy is conserved in each collision, then ln/(v) must 
be expressible as a linear combination of these collision invariants as 

ln/(v) = — (v - vo)^ + constant, (18) 

where G is the temperature in energy units {ks = 1) and m is the mass of a molecule. Here 
Vo is the velocity of the center of mass of the system. Hence we have shown that the system 
approaches the Maxwell-Boltzmann distribution. 

B. Structure of r(v^, v^; vcv), and connection with the Boltzmann equation 

We define new variables 

vt = (va + vb)/2, u = va-vb, m = |u| (19a) 
= (v + vc;)/2, u' = vc;-v, u'=|u'|, (19b) 



where vt and are the center of mass velocities before and after the colhsion. The Jacobian 
of the transformation is unity and integrations can be written in terms of the new variables. 
Momentum conservation is imposed on T(vyi, v^; vc, v) as 

r(vA, vb; vc, v) = 5%VT - v^)G(u, u'). (20) 

The integral in Eq. ( fTTl) is then written as 

I = j f{^rA)f{^B)T{wA, vs; Vc, v) d\A d\B dvc (21a) 
/(v + ^^)/(v - ^^)G(u, u') dudu'. (21b) 
The conditions in Eqs. ([1]) and ([2]) become for G{u,u'): 

j G{u, u')du = j G{u, u')du' = 1, (22) 
G(u,u') = ^(u^u). (23) 

Energy conservation requires that u = u' . If we define unit vectors u = u/u and n = u'/m 
and the angle 9 between them as (cos^ = u ■ fi), we can write G(u, u') as 

G(u,u') = ^^^^^^7(^,«). (24) 

The condition in Eq. ( l22l) becomes for g{6,u): 

jg{e,u)dh=l, (25) 

where g{d, u)dh is the probability of scattering into the solid angle dn in the center of mass 
frame. Then the integral / in Eq. (I21ap becomes 



/= Jf{^r+'^ + uh/2)f{v-'^ + uh/2)g{e,u)dudn. (26) 
If we write the /(v) term in Eq. (flTl) as 

/(v)= jf{^r)f{^r-u)g{9,u)dudh, (27) 
which follows from Eq. fl25|) and the normalization of /(v), we can write Eq. ffTTj) as 

^ J [/(v^)/(vb) - /(v)/(v - u)] g{e, u)dudh, (28) 



g/(v) 
dr 
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where 



u 

= v + - +Mn/2 
u 

Vfi = V - - + Mn/2, 



(29b) 



(29a) 



Equation (128|) is almost in the form of the Boltzmann equation. 

The Boltzmann equation represents a dilute gas for which the collision probability is 
proportional to uaxiu), where (7t{u) is the total cross section. We consider a large enough 
number R such that the ratio uaT{u)/R for a selected pair is almost always less than unity. 
For (Jt{u) = (To the constant R/ctq can be chosen to be a few (say five) times the rms velocity. 
Then when a pair is selected, we take a random number r and allow the collision to occur 
if r < uaT{u) / R\ we select another pair if r > uaT{u) / R. Although this procedure insures 
that the collision probability is proportional to U(Jt{u), it appears to violate the condition 
in Eq. fl25|) that all the selected pairs have a collision. To satisfy the condition in Eq. fl25|) 
we select g{6, u) as 



where a{6, u) is the differential cross section. The latter is related to the total cross section 
(Tt{u) by 



The second term in Eq. (!30l) transfers the initial velocities to the final velocities with the 
probability 1 — U(Jt{u)/R and the collision becomes a null collision. The (5(u — n) requires 
u = fi which implies u' = u since u' = u from the energy conservation and u' = u'fi, 
u =uvl. This means — = vc — v. We also have + = vc + v from center of mass 
velocity conservation. These two equations yield vc = and v = and velocities have 
not changed. A normal collision occurs with the probability U(Jt{u) / R. It is easy to verify 
that g{0,u) given in Eq. fl30|) satisfies the condition in Eq. fl25l) . 
If we substitute g{9,u) in Eq. (|30l) into Eq. (1281) . we obtain 



where and were given in Eq. (l29!l . Equation (132|) is essentially the Boltzmann equa- 
tion with the difference that /(v) is the probability density in velocity space whereas the 




(30) 




(31) 




(32) 
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Boltzmann equation is written in terms of the probability density in both physical and ve- 
locity space. If the volume of the cell containing the molecules is V, then we can write 
Eq. ([32]) for F(v) = (iV/1/)/(v) as 

= j [F(v^)F(vb) - F(v)F(v - u)] ua{9, u) dn dh, (33) 

where t = tV/RN = 2nV/RN'^ is interpreted as the physical time. Equation (l33l) is the 
Boltzmann equation for a homogeneous gas. 

III. DISCUSSION 

Let us summarize the direct simulation Monte Carlo algorithm for solving the Boltzmann 
equation. We choose a sufficiently large R such that only a negligible fraction of the selected 
pairs (say less than one in a thousand) violate the condition uaT{u)/R < 1. Then we select 
pairs randomly and let them collide with probability uaT{u)/R. The latter is achieved by 
generating a random number r and letting the collision occur if r < U(Jt{u)/R. If a pair 
collides, then in the center of mass system the collision occurs within the solid angle dn 
with probability P{6)dn = [o'{6,u)/o'T{u)]dh. Suppose that we put the z-axis along u and 
we need to determine n = u'/m, which is determined by the angles 6 and 0. To determine 6 
we need to generate a random value of 6 by converting the random numbers produced by a 
uniform probability distribution to random numbers in the interval (0, tt) according to the 
probability distribution P{6). The angles in the interval (0,2-7?) are equally likely. In this 
way we determine the final velocities of the particles as u' and — u' in the center of mass 
frame. By adding the center of mass velocity we find the final velocities in the lab frame. 
After storing the final velocities of the particles, we choose another pair and repeat the same 
process. The physical time is t = 2nV/N'^R, where n is the number pairs chosen to make 
attempts for a collision. If the number of collisions in a given time is required, we can count 



the successful attempts for a collision. In Ref. ^ this algorithm for keeping track of the time 
is called the "no time counter method." 

The original method of Bird^ to keep track of the time was the time counter method. 
Consider a narrow interval of uaxiu) values. For a large n there will be An pairs with 
U(Jt{u) values in this interval. Of these, only {u(Jt{u) / R)/S.n of them will make collisions 
corresponding to a time interval {2V/N'^R)An. Thus the elapsed time per successful attempt 
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IS 



^ {2V/N^R)An ^ 2V 

{uaT{u)/R)An N^uariu)' ^ ' 

In the time counter method we let every pair colhde, increase time hj At {t ^ t + At) 

after each colhsion, and keep selecting pairs and colliding them until we reach the desired 

time. Every collision will cause a different time increment depending on the value of uaxiu). 

One disadvantage of this method is that if a collision with a low U(Jt{u) occurs, the time 

increment will be large. Such collisions can occur with pairs having almost equal velocities. 

The time counter method was declared "obsolete' in Ref. sj. But it is useful to be aware of 

the method since it is widely used in the past and it might come across in some papers. 

If the purpose of the simulation is to demonstrate that the velocity distribution ap- 
proaches the Maxwell-Boltzmann distribution, we could let all the selected pairs make a 
collision and the velocity distribution will converge to the Maxwell-Boltzmann distribution. 
This simplification corresponds to u(7t{u)/R = 1 or axiu) = R/u, where the total cross 
section is inversely proportional to the relative velocity. Also, if it is desired to not discuss 
cross sections and the time tracking method, it is convenient to assume isotropic scattering 
in the center of mass frame. Then u' can be calculated by taking a random unit vector fi and 
multiplying it by m. These two simplifications make the programming easier and an under- 
graduate student with some programming background can write a program demonstrating 
the Maxwell-Boltzmann distribution. 

Direct simulation methods are also applicable to radiative processes and chemical reac- 
tions and the present formalism generalizes to all these cases in a more or less straightforward 
fashion for homogeneous gases. Such generalizations can be a useful teaching tool and a fer- 
tile field for student projects. 
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